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Abstract 

The effects of compressibility on Rayleigh-Taylor instability (RTI) are investigated by inspecting 
the interplay between thermodynamic and hydrodynamic non-equilibrium phenomena (TNE, HNE, 
respectively) via a discrete Boltzmann model (DBM). Two effective approaches are presented, one 
tracking the evolution of the local TNE effects and the other focussing on the evolution of the mean 
temperature of the fluid, to track the complex interfaces separating the bubble and the spike regions 
of the flow. It is found that, both the compressibility effects and the global TNE intensity show 
opposite trends in the initial and the later stages of the RTI. Compressibility delays the initial 
stage of RTI and accelerates the later stage. Meanwhile, the TNE characteristics are generally 
enhanced by the compressibility, especially in the later stage. The global or mean thermodynamic 
non-equilibrium indicators provide physical criteria to discriminate between the two stages of the 
RTI. 

PACS numbers: 47.ll.-j, 47.40.-x, 47.55.-t, 05.20.Dd 
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I. INTRODUCTION 


Rayleigh-Taylor instability (RTI) occurs at the interface between two fluids with different 
densities, subjected to an acceleration directed from the lower density fluid to the higher 
density one. A typical case is a heavy fluid resting on the top of a lighter one in the presence 
of a gravitational field. Under such conditions, density perturbations at the interface grow in 
time under the effect of gravity. The first detailed study of this instability was conducted by 


Rayleigh jl| in the 1880s. Later the first study was extended to accelerated fluids by Taylor 
^ in 1950. The first experiment was performed by Lewis in the evolution of an unsta 
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air-water interface [3|. Another experiment by Emmons et al. confirmed these findings 




Such an instability plays a prominent role in many natural and industrial processes, such as 


a, 


r. 


(ICE) i , type-la supernovae rJ, hot-wire diagnostics jSl], quantum magnetized plasmas 


colloidal mixtures 


3upei 

m 


devices__Jor sustainable energyj)roduction, say turbine_s 0 and inertial-confinement fusion 


0 , 


etc. 


In the above-mentioned fields, the compressibility effects on RTI are essential, and even 


dominate 


lll-ll4j|. deserving careful investigation. In fact, many t 


leoretical and numerical 


studies have been performed, especially on the initial linear stage 


15l-l2l|. In those studies, 


the compressibility effects on RTI growth rate are generally probed via changing the ratio 
of specific heats and the equilibrium pressure at the interface. Specifically, in 2007, Lafay 
et al. found that, in isothermal case, the stratification has a stabilizing effect while the 
compressibility has a destabilizing effect for two miscible viscous and compressible fluids 
j^. In 2008, He et al. reported that, in an inviscid case, the influences of the ratio 
of specific heats are as below. It mitigates the RTI when the upper heavy fluid is more 


compressible, while it enhances the RTI when the lower fluid is more compressible 


In 


2010, Ye et al. demonstrated that the compressibility has destabilizing effects for inviscid 
compressible fluid with exponentially variable density profile 0 

Although the compressibility effects have been studied extensively, several fundamental 


problems remain open, sue 
increasing compressibility 


1 as t 


le nonequilibrium effects in RTI, especially for the case of 
2ll-l28|. Eor the case with strong compressibility, the interfacial 
dynamics becomes more complicated as the RTI unfolds, resulting in very substantial gra¬ 
dient forces (Vp,Vu and VT) around the interfaces and very pronounced thermodynamic 
non-equilibrium (TNE) effects, where p,u and T are the local density, flow velocity and 
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temperature, respectively. The more pronounced the compressibility, the more complex the 
interfaces and the TNE effects as well. It is known that the Navier-Stokes model falls short 
of describing the complex interfaces and TNE effects At the same time, molecular 

dynamics and Monte Carlo simulations can not access macroscopic spatial-temporal scales of 
interest at affordable computational cost. Under such conditions, a kinetic approach based 
on a suitably simplified model Boltzmann equation is preferable. 



LB method was developed to probe the trans- and snpercritical fluid behaviors or both the 
hydrodynamic non-equilibrium (HNE) and TNE behaviours, which has brought some new 
physical insights into the fundamental mechanisms of the system. Physically, such an ex¬ 
tended LB kinetic model or discrete Boltzmann model (DBM) is roughly equivalent to a 
hydrodynamic model supplemented by a coarse-grained model of the TNE behaviours |59l |. 
The DBM has been applied in the combustion system, phase separation system, and com¬ 
pressible flow system with shocks j25l-l3l|. bnt still not in the RTI system. In this work, 
we further extend the DBM to investigate both the HNE and TNE behaviours in the com¬ 
pressible RTI system. Compared with previous studies on RTI, besides the compressibility 
effects, the interplays of varions non-eqnilibrium behavionrs are onr main concerns. 

The rest of the paper is structured as follows. In Sec. II, we first briefly review the DBM 
used in this work, then show the basic evolutions of so called “single” mode RTI and its 
TNE effects. A new front-tracking scheme based on the TNE property is presented in the 
same section. The effects of compressibility on RTI are studied in detail in Sec. III. Sec. IV 
snmmarizes and concludes the present paper. 
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II. EVOLUTIONS OF RTI AND ITS TNE CHARACTERIZATIONS 


A. Discrete Boltzmann Model 


Instead of using the traditional Navier-Stokes equations, in this work the compressible 
RTI system is described by the following discrete Boltzmann equation with Bhatnagar- 


Gross-Krook model 


m 




dfi ^ df. 


a • (v^ u) __ — (f. 

rji Ji \Ji Ji )) 

1 T 


( 1 ) 


where /j(r, Vj,t) is the discrete distribution function, r the spatial coordinate, t the time, Vj 
the discrete velocity model and i = 1,2^ ■ • ■ , N the direction of the discrete velocity model, 
u is the macroscopic velocity, a an external body force, r the relaxation time, and is the 
equilibrium distribution function. 


Following the ideas presented in Refs. 


25l-l3ll|. we use 


A* = M* (fA-M* (f'A 


( 2 ) 


to describe the TNE effects, where ^ represent the kinetic central moments, in which the 
variable v in Eqs. (IA4I) - (IA7I) (see the Appendix) is replaced by v* = v — u. M^, Mg, Mg 
and M 4 2 are associated with the Non-Organized Momentum Elux (NOME), Non-Organized 
Stress Flux (NOSF), Non-Organized Energy Flux (NOEF), and Flux of NOEF, respectively. 
Here, high order kinetic moments reflect the molecular individualism on top of organised 
collective motion, which we conventionally label as non-organised (NO) modes. 

In our simulations, we study the spatiotemporal evolutions of a single-component fluid 
when initially prepared on the hydrostatic unstable equilibrium, i.e., with a cold uniform 
region in the top half and a hot uniform region in the bottom half. In the two half volumes, 
we fix two different homogeneous temperatures, with the corresponding hydrostatic density 
profiles j^. We consider the lower and upper borders are kept far from the perturbation 
interface during the process of RTI, so there is no heat change in the lower and upper borders. 
Thus, adiabatic and non-slip boundary conditions are applied at the top and bottom walls 
and periodic boundary conditions on the horizontal boundaries. Specifically, at the top and 
bottom boundaries, we set the velocity to be zero. The forward Euler finite difference scheme 

n 

and the nonoscillatory nonfree dissipative scheme are used to discretize the temporal 
and spatial derivatives, respectively. 
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B. Evolutions of RTI 


The starting configuration of RTI is a compressible fiow in a 2D domain [—d/2, (i/2] x 
[—2d, 2d]. We consider two layers of the fluids at rest in the constant gravity field with 
the initial position of the interface ydx) = yocos{kx), where yo = 0.05d, the wave number 
k = 27:/X, and the A = d is the wavelength of the perturbation. The temperatures of two 
half volumes are initially constants and each of the half part is in hydrostatic equilibrium 

dyPoiy) = -gpoiy). ( 3 ) 

So the initial hydrostatic unstable configuration is given by 

Toiy) = Tu,poiy) = ^exp {‘2d - y)], y > y^x), 

< To{y) = Td,po{y) = ^exp [;^(2d - (4) 

d J-u 

-^{y - yc{x))],y < yc{x), 

K j-d 

where Pq is the initial pressure at the top boundary. and represent the initial temper¬ 
ature of the upper half part and the lower half part, respectively. Under this condition, we 
have the same pressure at the interface 


Pu^u Pd^d-) 


( 5 ) 


where and pd are the densities of the upper and the lower grid aside the interface. Then 
the initial Atwood number can be defined as Q 


^ ^ Tj T„ 

Pu X- Pd TdX- Tu 

Here we study both the hydrodynamic and thermodynamic behaviours of the single com¬ 
ponent compressible RTI. In our simulations, a grid size of 256 x 1024 is adopted. The other 
parameters are n = 3, c = 1.3, tjq = 15, Ax = Ay = 0.001, At = 2 x 10“^, Pq = 1.0, 
ttx = 0.0, tty = —g = —1.0, T = 5 X 10“^, Tu = 1.0 and At = 0.6. Figure H] displays the 
density evolutions of RTI. We observe that, at the beginning, thermal diffusion smoothes the 
discontinuous initial density interface, then a transition layer with finite thickness appears 
and the local effective Atwood number decreases. The amplitude of perturbation exponen¬ 
tially grows and the initial configuration remains cosine type until t = 0.6. After that, the 
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FIG. 1: (Color online) Density evolutions in the RTI simulated by the DBM at various times. 
The same results can also be obtained from the Navier-Stokes model. The larger the density, the 
stronger the inertial effects, the more difficult to change the velocity. Consequently, the upward 
perturbation grows into a shape similar to bubble, while the downward perturbation grows into a 
shape similar to spike. 

RTI enters the nonlinear stage, highlighted by the ontstanding spike and the appearance of 
Kelvin-Helmholtz instability due to the difference of the tangential velocity at the interface. 

Owing to the effects of gravity, the hngers of lighter fluid continuously penetrate into 
the heavier fluid, while the heavier fluid falls into the lighter one with a rolling-up process, 
resulted in the increase in the mixing layer amplitude, and forming a pair of secondary 
vortices which appear at the tails of the roll-ups, just like “mushroom” shape. The bubble 
rises due to the release of the compressive energy from the lighter fluid to the heavier 
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one. In the later stage, since the effects of the viscosity and thermal diffusion, the tails 
on both vortices gradually become less sharp and lon^narrow. The simulation results are 
qualitatively consistent with that of the experiments p, , reflecting the basic characteristics 
of the real physical process. 



FIG. 2: (Color online) Temperature profiles averaged in the x-direction versus y coordinate at 
different times. The mixing layer becomes wider with time. The perturbation depth in the upper 
part is smaller than that in the lower part. 

To quantitatively describe the characteristics of the mixing layer, we plot the averaged 
temperature profile T{y) against the y axis at t = 0.0, 0.2, 0.6, 1.0, 1.3, 1.6, 1.9, and 2.2 in 
Fig. [21 T{y) is defined as 

T{y) = j J^T{x,y)dx. (7) 

It varies from being discontinuous to being irregular that shows the thickness of the mixing 
layer and the amplitude of the temperature oscillation increases with time. The zig-zags in 
the profiles indicate the heat conduction of fluids from the high temperature area to the low 
temperature region and the irregularity in the mixing layer. 

C. TNE characterizations of RTI and corresponding interface-tracking technique 

Through the DBM, we can study not only the HNE behaviours, but also the TNE effects 
of RTI. The TNE effects can be interpreted as the manifestations of molecular thermo- 










fluctuations relative to the macroscopic flow velocity u and can therefore help in gaining a 
better understanding of the kinetic effects on the onset and development of the RTI. 
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FIG. 3: (Color online) Profiles of macroscopic qnantities and non-eqnilibrium effects along the 
central line x = Nx/2 at t = 1.6. The gradients of macroscopic qnantities work as driving forces 
of the TNE effects. The TNE quantities show specific TNE statns of the system. The most 
prononnced TNE effects occnr around the interface where the macroscopic quantities have large 
gradients. Because the gravity is in the y-direction, the “flux” in the y-direction is more pronounced. 
The non-zero valne of ^2xy ^ typical TNE qnantity. 

In Fig. 131 we illustrate the profiles of macroscopic and TNE quantities along the central 
line X = Nx/2 att= 1.6. From Fig. 131(a) one can observe that all the macroscopic quantities 
show complex behaviours around the interface. The temperature and density profiles show 
the largest and second largest gradients at this time, respectively. Non-equilibrium effects 
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are mostly pronounced around the interface, where the gradients of macroscopic quantities 
are particularly strong. Figure |3] (b) shows that the internal kinetic energies in the x and 
y degrees of freedom deviate from their equilibrium value oppositely in the same amplitude 
around the interface. Each one of and ^2yy deviates from zero oppositely in front of and 
behind the interface. /^2xyi which is zero at equilibrium, shows small but finite values around 
the interface, which is a typical TNE effect. From Figs. [3] (c)-(e) we can also appreciate 
that ^lyyy, ^^xxy ^ pBuks at the interface. Around the interface, 

^lyyy > 0, ^l^xy ^ O' ^hc valuc of /S.lyyy-\- iu Fig. |3] (c) can be read from the curve of 
A *3 ^^y in Fig. [HI (d). The positive peak of A *3 ^^y signals an upward internal kinetic energy 
flux. This is plausible, because heat transfers from higher to the lower temperature regions. 
In Fig. IHI (e), besides the non-zero value of ^*i^i2)xyi ^\i2)xx shows a larger amplitude than 

^(4,2)ro 

Usually, the depth of the mixing layer is an important parameter to measure the evolution 
of RTF We use the half amplitude to measure the mixing layer by capturing the spike and 
bubble. For incompressible RTI, this measurement is readily performed by tracing the 
constant density. However, in the compressible case, how to measure the mixing layer 
remains a thorny problem. Here we present two independent interface-tracking methods: 
(i) tracking the mean temperature of the upper and lower fluids, (ii) tracking the maximum 
values of TNE characteristic quantities, such as A *3 The second method is based on the 
fact that ^(3,1)2/ takes its maximum value at the position of interface along the ^-direction of 
the spike and bubble. We can adopt this TNE observable to capture the spike and bubble 
and obtain the thickness of the mixing layer, see Fig. jH The agreement between the results 
obtained from the above two approaches shows the effectiveness of the two tracking schemes, 
see Fig. [5l From Fig. |6l we find that the flow field is qualitatively consistent: from t = 0.2 
to t = 0.6, the perturbation amplitude grows exponentially, with a linear growth rate of 
about 0.082. 

III. EFFECTS OF COMPRESSIBILITY ON RTI 

According to the inviscid, isentropic Euler equation 

dtVL + u • Vu -F -c^Vp = g, (8) 

P 


10 



FIG. 4: (Color online) Positions of the bubble and spike versus time. The velocity of the bubble 
is smaller than that of the spike. The lighter fluid is “softer” and thus it is easier for the spike to 
pass through and grow. 



FIG. 5: (Color online) The perturbation amplitudes obtained by two different tracking approaches. 
The local TNE stength can be used to track interfaces. The good agreement shows that the two 
approaches validate each other. 


we can introduce a non-dimensional number Hi as below: 


lgl 

9 _ ^ 

YPc2 

P " 

kc^ k/g 


(9) 


It is clear that Hi can be regarded as the strength of the gravity relative to the gradient 
of pressure. Since dp/dp = describes the compressibility of the flow system, the non- 
dimensional parameter Hi can also be regarded as a relative compressibility. Since the 
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FIG. 6: (Color online) The growth of the perturbation amplitude by the DBM. The good linear 
fit is consistent with the linear theory on the RTF The result shows that the linear theory for 
incompressible flows works also for a period in the compressible flows. The DBM and the linear 
theory validate each other. 


compressibility is dimensional, it is not suitable for studying the effects of compressibility 
on RTI. Besides c*, both gravity g and the wave number k of the perturbation can influence 
the increasing rate of RTI. Under such conditions, Hi is a good non-dimensional parameter 
to describe the relative compressibility and is also a good parameter for studying the effects 
of compressibility on RTI. Hi is a variation of speed of sound and consequently a variation 
of 7 , and also of the stratification width, l/k. It can be controlled by adjusting g. However, 
it is totally unrelated to the viscosity and heat conduction. 

Similarly, 


H2 = T^/^ 


T 



( 10 ) 


can be regarded as the ratio of two time scales. Since the relation time r is relevant to the 
viscosity and thermal conductivity, the non-dimensional parameter H 2 can also be regarded 
as a relative viscosity or thermal conductivity. 

Therefore, we can define Hi = g/{kc^) and H 2 = to nondimensionalize the 

compressibility and the viscosity effects. The dimensionless time is then defined as t* = 
t/■\/27 t/{ kg) and the dimensionless length scale is 27r/fc. To study the compressibility Hi, 
we should make sure the viscosity H 2 is constant. Furthermore, we can set the initial 
T = 2.0 X 10“"^ and g = 1.0 to obtain the prescribed viscosity H 2 . So we can prescribe Hi by 
changing g. To ensure H 2 is a constant in all simulations, we need to recalculate r as derived 
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from H 2 = T^/^ each time. In our numerical simulations, we set the initial Atwood number 
At = 0.6, Tu = 1.0, Po = 1-0 uniformly for simplicity. Meanwhile, the mesh is specified by 
setting Ax = = 0.001, and Nx x Ny = 256 x 1024. Time step is At = 2 x 10“^. The 

other parameters used here are uniformly n = 3, c = 1.3 and rjo = 15.0. 



FIG. 7: (Color online) The growth of the dimensionless amplitude A* with different values of 
compressibility Hi. In the initial stage, the compressibility tends to inhibit the RTI, while in the 
later stage it tends to strengthen the RTI and the strengthening effect approaches saturation with 
increasing compressibility. 




FIG. 8: (Color online) The dimensionless amplitude versus the compressibility at (a) = 0.6 and 

(b) = 2.5. This figure shows more neatly the specific compressibility effects at one given time in 

each of the two stages. 

Figure [3 displays the time evolutions of the dimensionless amplitude A* with different 
values of the compressibility. It is found that the effects of compressibility can be roughly 
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divided into two stages: (i) In the first stage, for t* < 1.1, compressibility stabilizes the RTI; 
(ii) In the later stage, for t* > 1.1, compressibility accelerates the RTI. Particularly, we show 
the results at two dimensionless time instants t* = 0.6 and t* = 2.5, see Fig. |Hl which are 
fitted by two typical power-law relationships. When the compressibility is sufficiently large, 
its effects on the growth of RTI become much less evident. 

To interpret the above phenomena, we focus on the evolutions of the rates of the system 
internal and kinetic energies. For the rate of the internal energy, from Eq. flA8H (see the 
Appendix), we have 

P^ = -PV-u + V-(kVT) + Vu:P', (11) 

where d/dt = d/dt + u ■ V. The region that we measure lies in the middle of the system and 
its width is half of the system height. This region is large enough to prevent interference of 
the upper and lower boundaries with the spike and the bubble, i.e. no heat is supplied to or 
removed from the measured region as a result of any interaction with the boundaries. The 
heat conduction within the measured region makes no contribution to the increasing rate of 
internal energy. Therefore, when considering the rate of internal energy, the second term in 
the right-hand side of Eq. OB can safely be ignored. 

We define the rate of the compressive energy, Ec, as 

Ec = - [ P^- ndV, (12) 


where dV is the volume element, and the rate of the internal energy due to dissipation or 
viscosity E^ as 


E,= / Vu : P dV. 


(13) 


Eor the rate of the kinetic energy, from Eq. 

du 


(see the Appendix), we have 


p— = pa — VP -I- V ■ P . 
dt 


(14) 


This term includes three contributions. Similarly, we define the rate of the kinetic energy 
change by gravity, Ekg, as 

Ekg = y U • {pa)dV, (15) 

the rate of kinetic energy change due to dissipation, E^d, as 

Ekd = j u-(V-P')dR, (16) 
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and the rate of kinetic energy change by pressure, as 

Ekp = - ju-VPdV. (17) 

Figure El illustrates the time evolution of the various energy rates, Ec, E^, Ekg, E^d and 
Ekp, with different compressibilities. We first discuss observations for the second stage of 
RTI evolution. From Fig. IHl (a) it is clear that with decreasing the compressibility, Ec goes 
gradually back to the case of incompressible flows, where Ec = The extent of compressive 
energy rate Ec increases with time and compressibility. The rate Ec is negative, which means 
that the fluid volume expands in time. Since the fluid has internal dissipation, the rate of 
the energy dissipation Ed increases in time and with increasing viscosity [see Fig. El (b)]. 
In this study, the compressibility Hi is increased by increasing the gravity acceleration g, 
consequently the kinetic energy rate by gravity increases with increasing compressibility Hi 
[see Fig. [9] (c)]. Since the volume expansion rate increases with increasing compressibility 
and time, the rate of kinetic energy by dissipation or viscosity Ekd show similar behaviour 
[see Fig. El (d)]. It is interesting to observe that the rate of kinetic energy by pressure E^p also 
increases with the compressibility Hi [see Fig. El(e)]. This behaviour can be understood by 
considering that the compressibility Hi increase means that gravity acceleration g increases. 
Consequently, the pressure gradient, VP, becomes larger, thereby enhancing the rate of 
kinetic energy by pressure Ekp as well [see Eq. fTTZl) ]. 

Now, we come back to the initial stage. Initially, the flow velocity is low, and so is the 
Mach number, so that the compressibility of the fluid is naturally small. So, the values 
of Ec, Ekg and Ekp are also small. The initial state of the system around the interface is 
thermodynamically unstable. Without any perturbation, the interface molecules are in a 
mechanically stable but thermodynamically unstable state. In this case, the forces due to 
the temperature and density gradients tend to increase the boundary layer and decrease 
the local Atwood number. If the interface is initially perturbed, most interface molecules, 
in between the crest and trough, experience a gradient force with a non-zero horizontal 
component, which tends to flatten the interface. It should be pointed out that in our 
numerical simulations, the initial state is slightly different from the one given by Eq. (|H), 
due to the finite lattice spacing and the free-energy of the numerical initial state is higher 
than the one given by Eq. (|H). The system relaxes towards its minimum free energy state 
and gravitational energy transforms into kinetic energy. Due to viscosity, while approaching 
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FIG. 9: (Color online) The time evolutions of the various energy rates at different values of com¬ 
pressibility Hi. The energy rates include the rates of compressive energy Ec^ internal dissipation 
energy kinetic energy induced by gravity Ekg^ kinetic energy induced by dissipation Ej^g and 
kinetic energy induced by pressure E]^p. It is clear that a larger compressive energy rate is involved 
in the later stage. This figure highlights significant differences as compared to the incompressible 
scenario. 

this minimum free energy state, the local flow velocity tends to zero, so that the amplitudes 
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of Ed and E^di after a quick initial rise, decrease significantly in the longer term. The stage 
for the initial quick increase is indeed very short. 

The process of RTI can also be divided into two stages and be interpreted from the 
point view of energy transformation: (i) In the initial stage, the compressibility provides 
a stabilizing effect. This is mainly because a higher compressibility Hi corresponds to a 
larger gravity acceleration corresponding in turn to a larger local density p or pressure 
P, and consequently a higher heat conductivity k. The heat conduction tends to decrease 
the local Atwood number and broaden the interfaces of the density and temperature, (ii) In 
the later stage, the compressibility has a destabilizing effect, which can be explained as the 
transformation of the stored compressive energy into kinetic energy. 

To provide an estimate of the TNE effects resulting from compressibility, we follow the 
same idea used in Refs. Q and j^, and define a global average non-equilibrium intensity 
or “TNE strength” 

D* = + (A^)2 + (A* (18) 

where A)„^^ is the average absolute value of TNE components. 



FIG. 10: (Color online) The time evolntion of the global average TNE strength with different valnes 
of compressibility Hi. The compressibility first decreases in the first stage and then increases in 
the later stage the global average TNE strength. The higher the compressibility, the stronger the 
global average TNE effects. 

Figure [1(11 shows how the compressibility affects the global average D*. With increasing 
compressibility, the global average deviation from thermodynamic equilibrium increases. 


17 





Since the initial condition is in thermal non-equilibrium state, the system has a tendency 
to approach the thermodynamic equilibrium state at hrst. Therefore, in the first stage, 
D* shows a decreasing trend. In the later stage, D* is found to grow exponentially. This is 
because TNE effects are tightly coupled to the interface dynamics, which shows an increasing 
area (coarsening) and morphological complexity. 



FIG. 11: (Color online) The average valnes of various components of TNE quantities versus com¬ 
pressibility at four nondimensional times: t* = 0.1, 0.6, 2.0 and 2.5. The fignre shows more 
specific global or mean TNE effects than Eig. [101 The relative strengths of those dependences on 
compressibility may change with time. 


Figure rm shows the average value of each component of TNE observables with changing 
compressibility at four different times t* = 0.1, 0.6, 2.0 and 2.5. It is found that: (i) All 
TNE dynamic modes increase as the compressibility increases, at all times. Since all the 
observations have a limited accuracy, a dynamic mode is not visible whenever its amplitude 
or strength is smaller than a critical value, say 10“^. Therefore, one can observe that 
more TNE dynamic modes emerge and stand out with increasing the compressibility. At 
later stage, the higher order terms of the TNE dynamic modes, such as ^* 4 ^ 2 )xy' ^ 

more important role. A single TNE dynamic mode is not the only decisive factor for the 
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amplitude of the D*; (ii) The relative strength among the dynamic modes may change at 
different stages, for example, is less than at the first stage, but reverse at 

the later stage; (hi) The strengths of some dynamic modes are always relatively small, such 
as ^ 2 xy ^ 2 yy Thcsc details are complementary with the above highly synthetic 

“TNE strength” D*. 

IV. CONCLUSIONS 

The Rayleigh-Taylor instability in compressible flows is studied via a discrete Boltzmann 
model. Besides hydrodynamic behaviour, the thermodynamic non-equilibrium effects most 
relevant to the hydrodynamic behaviour have also been studies in much detail, up to Atwood 
numbers around 0.9. It is found that the process of the Rayleigh-Taylor instability in com¬ 
pressible flows can be divided into two stages, exhibiting opposite compressibility effects: in 
the initial stage, compressibility stabilizes the Rayleigh-Taylor instability, while in the later 
stage it accelerates it. The physical reasons are as follows. A higher compressibility leads to a 
stronger gravity acceleration, which corresponds to a higher local pressure and consequently 
to a higher heat conductivity. In the first stage, the heat conduction tends to decrease the 
local Atwood number and broaden the interfaces of the density and temperature profiles. In 
the second stage, part of the compressive energy stored in the fluid is released and trans¬ 
formed into kinetic energy, thereby accelerating the Rayleigh-Taylor instability. The local 
thermodynamic non-equilibrium indicators provide useful observables to physically track the 
interfaces. Indeed, the global or mean thermodynamic non-equilibrium indicators permit to 
discriminate the two stages of the Rayleigh-Taylor instability. In the first stage, the system 
slowly evolves towards its equilibrium, while in the later stage, as the interface develops, the 
system moves away from local equilibrium, especially in the regions near complex interfaces. 
The above behaviour is enhanced by increasing compressibility, and so are the amplitudes 
of thermodynamic non-equilibrium kinetic modes. Besides a deeper physical insight into the 
kinetic procedures, the methodology and resulting observations may help to formulate more 
accurate meso and macroscale models for the complex compressible phenomena. 
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Appendix A 


In this work we use a 2D 16-velocity (D2V16) model schematically shown in Fig. lAll 
(|vj| = c for i = 1, • • • , 4, |vj| = 2c for i = 5,..., 8, |vj| = 3c for 9, • • • , 12, |vj| = 4c for 
i = 13, • • • ,16. = r]o for z = 1, • • • ,4, and ? 7 j = 0 for i = 5, • • • , 16.) 



FIG. Al: Sketch of the DVM used in the present paper. The nnmbers in the fignre are the indexes 
of the discrete velocities. 


The equilibrium distribution function satisfy the following seven kinetic moments [56] 


Mo(/r) = = »• 

i 

M:(/'’) = E 

i 


(Al) 

(A2) 
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M 2 ,o(/r) = = p[{D + n)T + u-u], (A3) 

i 

M2(D = E= + (A'*) 

i 

M 3 ,i(/D = X] = pu[{D + n + 2)T + u ■ u], (A5) 

i 

M 3 (/®'^) = ^ /rviVjVj = p[r(u„e; 3 e^( 5 ^^ + u^e„e^(5„^ + n^Bc^epScp) + uuu], (A 6 ) 


M4,2(/r) = E/r(v.-v, + r^,V.v. 


(A7) 


— pT \_{D -\- n -\- 2')T + u • u] Gct^fjSap + puu \_{D + 77. + 4)T' + u • u], 

where D is the space dimension. Besides the translational degrees of freedom, ry is a free 
parameter introduced to describe the n extra degrees of freedom corresponding to the molec¬ 
ular rotation and/or vibration. stands for that the m-th order tensor is contracted 

to a n-th order one. For 2D case, Mq and are scalars; and are vectors; 
and M 4 2 are the 2nd order tensors; Mg is a 3rd order tensor. The trace of moment is 
a conserved quantity and describes the energy, and its off-diagonal components correspond 
to the shear effects, which may not be zero whenever the system is in non-thermodynamic 


equilibrium state 


28l |. From the Chapman-Enskog analysis, such a DBM can recover the 


following Navier-Stokes equations 

dtp + V ■ (pu) = 0, 


dt{pu) -F V • (PI + puu) - V • P' = pa. 


a 


p(e + -|u| 


+ V 


p(e+ -|up )u -I- Pu 


(AS) 


with 


= V • ( kVT -|- u • P' -|- pu • a, 


P' = p Vu -I- (Vu)^ — 


D + n 


V-uI 


(A9) 


in the continuum limit, where 


^ D + 

P = pT, e = -^—T 


(AlO) 
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stand for the pressure and the total internal energy, respectively, /i, k are the dynamic 
viscosity and thermal conductivity coefficients having the following relations with r 

/i = Pr, n = CpPr (All) 

where Cp is the specific-heat at constant pressure. 
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